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We demonstrate a novel algorithm for computing the sensitivity of statistics in chaotic 
flow simulations to parameter perturbations. The algorithm is non-intrusive but requires 
exposing an interface. Based on the principle of shadowing in dynamical systems, this 
algorithm is designed to reduce the effect of the sampling error in computing sensitivity 
of statistics in chaotic simulations. We compare the effectiveness of this method to that of 
the conventional finite difference method. 


Nomenclature 


m x 1 vectors of ones 
Speed of sound 
m x 1 vector of homogeneous tangent solution weights 
m x 1 vector with projection of tangent onto homogeneous tangent solutions 
Drag coefficient 
Lift coefficient 
Drag 
Design variable of interest 
Objective functions or instantaneous contribution to objective functions 
Objective function sensitivities 
n x 1 linearization of instantaneous contribution to objective function 
Time dilation instantaneous contribution to sensitivity of objective function 
Number of time segments 
Lift 
m x m identity matrix 
Mach number 
Number of homogeneous tangent solutions 
Number of states (degrees of freedom) 
nm xX 1 vector of conserved variables 
n xX 1 spatial residual vector 
8,b Lorenz system parameters 
Time 


ry 
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T Time-averaging window length 

U Fluid velocity 

Vv nx m tangent instantaneous solution matrix 
x,y,z Lorenz system states 

At Time step 

€ Finite difference perturbation 

p Fluid density 

® Map to state at time t 


WwW Objective-function windowing weight 
Subscript 

j jth column of a matrix 

lo) Freestream value 

Superscript 

P Final state or objective in a time segment 
I Initial state or objective in a time segment 
a Time segment 


Indicates time-averaged quantity 


I. Introduction 


Many important phenomena in aerospace engineering applications exhibit chaotic dynamics. Examples 
of chaotic simulations include high-fidelity, scale-resolving turbulence simulations,! such as direct numerical 
simulations (DNS), large eddy simulations (LES), or detached eddy simulations (DES). Aerospace applica- 
tions including jet engines,” scramjet engine combustors,? and complex launch vehicle configurations* need to 
be analyzed with DES or LES to capture key physical phenomena such as flow separation. Multidisciplinary 
simulations, such as coupled fluid-structure simulations, can also exhibit chaotic dynamics.” * 

Design optimization requires computing the sensitivity of quantities of interest with respect to design 
parameters. Such sensitivity analysis is particularly challenging when the quantity of interest is an infinite- 
time-average, also called a statistic, of quasi-equilibirium chaotic dynamics. One example of such a quantity 
of interest is the long-time-averaged root-mean-square (RMS) norm of pressure on a launch vehicle.4 A 
finite-time-average approximation of any statistic, necessarily incurs a sampling error,® 


eT = ol aus (1) 


where T is the time averaging window size and c is some constant. This sampling error decays with the 
square root of the averaging window and makes the conventional finite-difference method very noisy for most 
practical unsteady simulations. 

To illustrate the impact of the sampling error on sensitivities computed with finite-differences, we use 
the Lorenz attractor, a simple chaotic system with three degrees of freedom 


WY _ ap — 2) ty =0 (2) 
—-a2x2“(r-z = 

7 y =0, 

 — ay + be =0. 


The sensitivity of the time averaged objective 
1 TotT ‘ 
= (z — 28)* dt, 
T To 


with respect to the Rayleigh number r is computed. Figure 1 (a) plots the objective function, averaged 
over T = 5 time units following To = 5 time units of run up time. Note that the longest time scale of the 


2 of 8 


American Institute of Aeronautics and Astronautics 


Objective function 


110 


105} 


100} 


95} 


90} 


85} 


80} 


75} 


70} 


65 


27 


28 29 30 31 32 33 
Design parameter 


34 


35 


Design parameter 


(c) T = 500 
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(d) T = 5,000 
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Figure 1: Time average of (z — 28)? in the Lorenz attractor (Eq. (2)) over four different time spans T, as a 
function of r, for b = 8/3 and s = 10. 


dynamics of the Lorenz attractor is about 1 time unit. Despite being averaged over a time span several times 
longer than the longest time scale of the system, the computed statistic shows no trend. The noise due to 
sampling error is so large that it is impossible to infer how the objective function depends on the parameter. 

To reduce the sampling error, one can increase the averaging length, as shown in Figures 1 (b), (c), and 
(d). These figures indicate that there is a clear minimum of the statistic when r ~ 31. From Eq. (1), 
reducing sampling error by a factor of 10 would require increasing the window size by a factor of 100. For 
high-fidelity simulations, the computation time becomes infeasible fairly quickly. 

The least-squares shadowing (LSS) method® has been developed to address this challenge in sensitivity 
analysis. The original LSS algorithm was intrusive, requiring significant modifications to any nonlinear 
solver, including development of an appropriate strong linear solver. A less intrusive algorithm has recently 
been developed.!° This paper demonstrates its utility on a chaotic, eddy-resolving turbulent flow simulation. 


II. The Non-Intrusive Shadowing Algorithm 


The algorithm presented in this paper is a modification of the non-intrusive Least Squares Shadowing 
(NI-LSS) algorithm’? that aims to further reduce the intrusiveness. In particular, the algorithm presented 
in this paper requires neither a tangent (linearized) solver nor an adjoint solver, approximating the tangent 
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algorithm presented by Ni et al.1° with finite differences. 


We assume that we have a simulator that solves an n-dimensional autonomous ODE, 
dQ 


Or +R(Q,D) =0, (3) 


for a given initial state, Q’, a given parameter D, and a given time t. The simulator would compute the 
solution state after the ODE evolves for time t. We also assume that the simulator computes a number 
of quantities of interest, F(Q,D), averaged over the time span t. Lastly, we also need it to compute these 
quantities of interest at the end of the time integration. 

Mathematically, the simulator can be represented as three maps: 


1. 6(Q,D,t) gives the solution state after a time t evolution of the ODE at parameter D, starting from 
initial state Q. 


2. F(Q,D,t) gives the time averaged F(Q,D) over a time t evolution of the ODE at parameter D, 
starting from initial state Q. 


3. F*(Q,D,t) gives F(®(Q, D,t),D), the instantaneous value of the quantities of interest at the end of 
the evolution. 


The NI-LSS algorithm runs this simulation to perform sensitivity analysis of the long time averaged 
objective F, 
1 7? 
f:= lim z/ F(Q, D)dt, 
T Jo 


T—-00 
with respect to D. It does this by computing m + 1 finite-difference approximated tangent solutions in K 
size T/K time segments, then solving a size mn x mn least squares system. The algorithm proceeds as 
follows: 


1. Choose an integer m, an estimate of the number of unstable modes of the ODE. Choose a time span 
At, an estimate to the reciprocal of the largest Lyapunov exponent.!! Then choose the number of time 
spans T/At, such that T is a multiple of the largest time scale in the system. 


2. Run the simulation for sufficient time to a quasi-equilibrium-state. Denote the state as Q°. From a 
dynamical systems point of view, this means that Q®° can be considered as being on the attractor of 
the system. 


3. Set V° to a random matrix of dimension (n,m) where n is the number of degrees of freedom of the 
dynamical system. Set V° to a zero vector of dimension n. Start from time segment i = 1. 


4. Compute R(Q*!,D). This can be obtained either directly from the simulator, or approximating it 
with —(®(Q‘“!, D, 6t) — Qi) /dt » —% for a small dt. 


5. Project each column of V‘~! as well as V‘~! against R, ice., 


RRivi-! ae Lawes RRivi-1 


i-1 _ yi-1 
™ " R™R °° R?R ° 


6. Run the simulator to compute 


Qi =8(Q'.D,At), F =F(Q'?,D,At), F®'=F*(Q*),D, At), 


7. For j =1,...,m, denote V% as the jth column of the n x m matrix V’. Run the simulator to compute 
 &(Q1 4 eVE_D, At) — Q F(Q'"1 + Vi}, D, At) -F Rrvi 
Vi= (Q J ) : ’ G; = (Q : ’ Hi; = ’ 
: € € AtR?R 


8. Run the simulator to compute 


vig ®(u;-14+ Eve" De, At) — Q! Giz F(Qu t+ eVi-1 De, At) — F ie R'vi-} 
a € ’ 7 € : * AtRTR’ 
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9. Compute the QR decomposition 
vi = OR’, 


then compute 
Bi os orvi 


and set V; = Q, V' = V' — QB’. 
10. Loop back to Step 4 with i =i+1, untili > T/At= K. 


11. Solve the reduced Least Squares Shadowing system 


x -R} I | - | B! 
min - s.t. feo. gh - as Wea 
A A 
-R* I B* 
A +1 5 A*kt1 
compute the sensitivity 
df 


K 
© wi ((Gi + Hat + (G' +4), 
i=l 


where w; is an averaging weight series satisfying yaa w' = 1. In this paper, 


ae 
ae (sin ea) 

Kk+1 

This algorithm is implemented in the open source repository https://github.com/qiqi/fds . All the results 
presented in this paper are based on this implementation of the algorithm. 


III. Two-Dimensional Chaotic Airfoil 


To demonstrate the NI-LSS algorithm, we consider a NACA 0012 airfoil at an angle-of-attack of a = 20° 
with a freestream Mach number M,, = 0.1 and a Reynolds number of 10,000. This flow is similar to the 
chaotic vortex shedding flows previously studied by Pulliam.'? The flow was simulated using FUN3D’s 
compressible, second order, finite volume flow solver. 


A. Sensitivity Estimates 


To demonstrate the NI-LSS algorithm, we consider the sensitivity of time-averaged lift and drag to pertur- 
bations of the freestream Mach number, M., and the angle of attack, a. The Mach number sensitivities can 
be estimated analytically by assuming that dC, /dM. = dCp/dM, = 0. This is a reasonable assumption 
for the small Mach number, M,, = 0.1, considered in this paper. Since the flow is nearly incompressible the 
effect of Mach number on Cz, and Cp is O(M7?).14 

If dC, /dM,. = 0, then for a unit reference area 


1 1 
L= 5 PoUocCL a 5 Poo too Miso 


Since the freestream density px. and speed of sound ag are constants 


dL > L 
3 pa MO = 4 
dM,,  °~%~ CL M.. (4) 
By the same reasoning, 
dD 2 D 
— = 0~05,MaCp = 2—— 5 
Nig PE °) 


Figure 2, shows that the analytical estimates are quite accurate. 
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Figure 2: Time-averaged lift and drag versus Mach number. The dotted line is a linear regression, the solid 
lines show local gradients computed with Eqs. (4) and (5). All cases were run for 500,000 time steps and 
averaged for 490,000 steps. Error bars are estimated using the batch mean approach with four means.!% 
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Figure 3: Time-averaged lift and drag versus angle of attack. All cases were run for 1,010,000 time steps and 
averaged for 1,000,000 steps. Error bars are estimated using the batch mean approach with four means.!* 


We do not have similar analytical estimates for the sensitivities to the angle of attack. In an attempt to 
numerically estimate the sensitivities, we performed 11 simulations in the range of a € [19.95, 20.05], each 
for 1,010,000 time steps. The lift and drag were averaged over the last 1,000,000 time steps, corresponding 
to 100,000 non-dimensional time units (non-dimensionalized with respect to the freestream speed of sound 
and the chord length), or 10,000 flow-through times (non-dimensionalized with respect to the freestream 
velocity and the chord length). The time-averaged lift and drag are shown in Figure 3. In spite of these 
lengthy calculations, no definite trend can be inferred. 


B. LSS Sensitivity Results 


The NI-LSS algorithm has been applied to estimate sensitivities with respect to the Mach number and the 
angle of attack. One hundred time segments are used in each calculation. Each time segment contains 200 
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Figure 4: Shadowing sensitivity for the derivative of lift and drag to the Mach number. The x-axis indicates 
the number of time segments, each containing 200 time steps, used in the shadowing sensitivity analysis. 
The light-colored bands indicate the analytical estimates. 


time steps, which corresponds to 20 non-dimensional time units (non-dimensionalized with respect to the 
freestream speed of sound and the chord length), or two flow-throughs (non-dimensionalized with respect 
to the freestream velocity and the chord length). In the calculation, we used m = 16 homogeneous tangent 
solutions, which requires m+ 2 = 18 simulations for each time segment. The total computation involves 
18 x 100 x 200 = 360, 000 time steps. The sensitivities are computed after each time segment is completed. 
Figure 4 shows the lift and drag sensitivities with respect to the Mach number versus the number of time 
segments. Figure 5 shows corresponding sensitivities with respect to the angle of attack. 

In Figure 4, we can see that the NI-LSS sensitivities approach the analytical estimates as the number of 
time segments increases. However, we do not observe a similar convergence for the sensitivity with respect 
to the angle of attack. We think that the lack of convergence arises because the sensitivities with respect to 
the angle of attack are small; and as a result, neither finite difference (as shown in Figure 3) nor the current 
NI-LSS algorithm (as shown in Figure 5) can accurately compute such small sensitivities. 


IV. Conclusion 


We presented a modified non-intrusive least-squares shadowing (NI-LSS) algorithm for computing the 
sensitivity of statistics in chaotic flow simulations. The algorithm minimizes modifications of the flow solver. 
The NI-LSS algorithm is applied to a two-dimensional simulation of a stalled NACA 0012 airfoil at Reynolds 
number 10,000, performed using the FUN3D flow solver. Sensitivities of drag and lift with respect to the 
Mach number and the angle of attack have been computed and compared to finite-difference estimates. 
Mach number sensitivities have also been compared with analytical estimates. The NI-LSS sensitivities 
with respect to Mach number agree with finite-difference and analytical estimates. Neither NI-LSS nor the 
conventional finite-difference approach is able to estimate the sensitivities with respect to the angle of attack. 
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Figure 5: Shadowing sensitivity for the derivative of lift and drag to the angle of attack. The x-axis indicates 
the number of time segments, each containing 200 time steps, used in the shadowing sensitivity analysis. 
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